plus sc10x

load dependancies

load 10x data

GEX.seur <- readRDS("./sc10x_LYN.marker.sort0829.rds")
GEX.seur
## An object of class Seurat 
## 17590 features across 24766 samples within 1 assay 
## Active assay: RNA (17590 features, 1199 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
color.FB <- ggsci::pal_igv("default")(49)[c(8,33,40,
                                            34,26,1,28,
                                            2,43,18)]
color.cnt <- color.FB[c(3,2,1,
                        10,9,8,
                        7,6,4)]
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "sort_clusters")

cleaning up

keep only Microglia clusters
remove DF0.05 doublets
remove AD.GFP.2

GEX.seur <- subset(GEX.seur, subset = preAnno %in% c("Microglia") & 
                                      DoubletFinder0.05=="Singlet" & 
                                      FB.info %in% setdiff(levels(FB.info),"AD.GFP.2"))
GEX.seur <- subset(GEX.seur, subset = percent.mt < 8 & nFeature_RNA < 2800 & nCount_RNA < 8500)
GEX.seur
## An object of class Seurat 
## 17590 features across 21021 samples within 1 assay 
## Active assay: RNA (17590 features, 1199 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.info", cols = color.FB[c(1,2,3,4,6,7,8,9,10)])

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "sort_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "sort_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "cnt", cols = color.cnt)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "cnt", cols = color.cnt)

plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "cnt", cols = color.cnt) 
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "cnt", cols = color.cnt) 
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "cnt", cols = color.cnt) 
plota + plotb + plotc

par(mar=c(6,4,2,3))
sl_stat <- table(GEX.seur$cnt)
barplot(sl_stat,ylim = c(0,3800),
        #col = c("#FF6C91","lightgrey",color.FB),
        col = c(color.cnt),
        main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:9*1.2-0.48 ,levels(GEX.seur@meta.data$cnt), las=3, cex.axis=0.85)
text(x=1:9*1.2-0.45,y=sl_stat+245,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)

new

re-clustering

GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1800)

# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2

head(VariableFeatures(GEX.seur), 300)
##   [1] "Cxcl10"        "Spp1"          "Ccl4"          "Cxcl13"       
##   [5] "Ccl5"          "Ccl3"          "Cxcl9"         "Ccl12"        
##   [9] "Spp1-EGFP"     "Mmp12"         "Cd74"          "H2-Aa"        
##  [13] "Gpnmb"         "Ccl2"          "Rsad2"         "Cst7"         
##  [17] "Cck"           "Top2a"         "H2-Ab1"        "Lpl"          
##  [21] "Ttr"           "Lgals3"        "Mki67"         "Fabp5"        
##  [25] "Iigp1"         "Ifi27l2a"      "Fn1"           "H2-Eb1"       
##  [29] "Apoc1"         "Ccl6"          "Gdf15"         "Saa3"         
##  [33] "Egr1"          "Ifitm3"        "Apoe"          "Enpp2"        
##  [37] "Wfdc17"        "Ifit3"         "Ch25h"         "Il1b"         
##  [41] "Prc1"          "Il1rn"         "Lyz2"          "Hist1h1b"     
##  [45] "Tnf"           "Ifit1"         "Ccl7"          "Stmn1"        
##  [49] "Pclaf"         "Cxcl2"         "Egr3"          "Ifit2"        
##  [53] "Igf1"          "Cenpf"         "Nusap1"        "Ube2c"        
##  [57] "Gm1673"        "Dqx1"          "Cd5l"          "Aldh1a3"      
##  [61] "Igfbp5"        "Ank"           "Cdkn1a"        "Ly6k"         
##  [65] "Apoc4"         "Gm43409"       "Isg15"         "Slc7a11"      
##  [69] "Wdcp"          "Ms4a7"         "Clec4d"        "Hmmr"         
##  [73] "Itgax"         "Prdx1"         "Cox6a2"        "Cd52"         
##  [77] "Ccl8"          "Csf1"          "Birc5"         "Atf3"         
##  [81] "Ctla2a"        "Ifnb1"         "Ctsd"          "Cd72"         
##  [85] "Syngr1"        "Cybb"          "Gbp2"          "Gnas"         
##  [89] "Vim"           "Ctsb"          "Bst2"          "Ly6a"         
##  [93] "Postn"         "Olfr1369-ps1"  "Lilrb4a"       "Tyrobp"       
##  [97] "B230312C02Rik" "Mif"           "Cd69"          "Neil3"        
## [101] "Gm45184"       "Aplp2"         "Plp1"          "C4b"          
## [105] "Dennd5b"       "Gm26885"       "Nr4a1"         "Fth1"         
## [109] "Edn1"          "Ccnb2"         "Lgals1"        "Ifit3b"       
## [113] "Ftl1"          "F830016B08Rik" "Cd63"          "Olr1"         
## [117] "Ramp1"         "Pbld1"         "Lgals3bp"      "Chac1"        
## [121] "Cdca8"         "Bex1"          "Pbk"           "Cenpa"        
## [125] "Slfn5"         "Trib3"         "Ifi204"        "Clu"          
## [129] "Hmgb2"         "Knl1"          "Hif1a"         "Ccr1"         
## [133] "Hist1h3c"      "Flt1"          "Ifi207"        "Serpine2"     
## [137] "Mamdc2"        "Gadd45b"       "Ctsl"          "Ifi44"        
## [141] "Pou3f1"        "Gria2"         "Sprr1a"        "Oasl1"        
## [145] "Ccl9"          "Colec12"       "Ctsz"          "Ifi211"       
## [149] "Tcap"          "Kif11"         "Id3"           "Fcgr4"        
## [153] "Ifitm1"        "Ccdc63"        "Ms4a4b"        "Id2"          
## [157] "Esco2"         "Dkk2"          "B930036N10Rik" "Il4i1"        
## [161] "Gdf3"          "Grb14"         "Gpr85"         "Rgs1"         
## [165] "Irf7"          "Tpx2"          "Rrm2"          "Ptgds"        
## [169] "Fxyd5"         "Enah"          "Cstb"          "Gm26737"      
## [173] "Lgi2"          "Lockd"         "Phlda3"        "Trem2"        
## [177] "Cxcl14"        "Cd9"           "Acod1"         "Baiap2l2"     
## [181] "Hist1h2ap"     "Gm34455"       "Crip1"         "Ecrg4"        
## [185] "Cenpe"         "Mt1"           "Oasl2"         "Timp2"        
## [189] "Ccl22"         "Clec12a"       "Egfl7"         "Rhox5"        
## [193] "Spag5"         "H2-Q7"         "Gapdh"         "Loxl2"        
## [197] "mt-Co1"        "Slamf7"        "Tspo"          "Gm21738"      
## [201] "Gm48837"       "Uchl1"         "Ccna2"         "Stap1"        
## [205] "Ptprf"         "Ppp1r14a"      "Alpk2"         "Crlf2"        
## [209] "Plac8"         "Aspm"          "Gm8251"        "Rab3c"        
## [213] "Ccnb1"         "Hbegf"         "Capg"          "Kif23"        
## [217] "Gfod2"         "Rpl41"         "Usp18"         "Snca"         
## [221] "Ifit1bl1"      "Sgo2a"         "Gm10800"       "Mgl2"         
## [225] "Shcbp1"        "Gpr84"         "Pkm"           "Hp"           
## [229] "2810403D21Rik" "Gm49756"       "Cdc20"         "Ifi208"       
## [233] "Kcnk9"         "Pld3"          "Fabp3"         "Plk2"         
## [237] "Cspg4"         "Nes"           "Spc24"         "Cd40"         
## [241] "Cd83"          "Kcnq1ot1"      "Gm10457"       "A330032B11Rik"
## [245] "Scd2"          "Ankle1"        "Gm4951"        "Slc7a14"      
## [249] "Tlr2"          "Anln"          "Tnip3"         "Rcan1"        
## [253] "Ms4a4c"        "Tk1"           "Slfn2"         "Ptgs2"        
## [257] "Gm17251"       "Tbx19"         "Meis2"         "Lrrc31"       
## [261] "Sh3d19"        "Rgs17"         "Olfr889"       "Stat1"        
## [265] "Peg3"          "Tnfrsf11b"     "Ptprz1"        "Axl"          
## [269] "Zbp1"          "Madcam1"       "Lhfpl4"        "Apod"         
## [273] "H2-D1"         "Adamts1"       "Mal"           "Plau"         
## [277] "Pf4"           "Atp6v0d2"      "Cd68"          "Syt5"         
## [281] "Selenow"       "Gm4841"        "H2-K1"         "1700120C14Rik"
## [285] "Gm29773"       "Slc1a2"        "Gm14636"       "Fam20c"       
## [289] "Serpina3g"     "Cd36"          "Bub1"          "Cit"          
## [293] "Iqgap1"        "Aurka"         "Smc2"          "Pcgf2"        
## [297] "Ptgs2os2"      "Smtnl1"        "Mir670hg"      "Gchfr"
# exclude MT genes  and more 
#   add sex-related Xist/Tsix
DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|Egr1|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AC|^AI|^AA|^AW|^AY|^BC|^Gm|^Hist|Rik$|-ps|Xist|Tsix|^Ifi|^Isg|^Mcm",
            rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
MT_gene <- grep("^mt-",rownames(GEX.seur),value = T)

VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
                                                c(MT_gene,
                                                  DIG, 
                                                  CC_gene) )

GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur), seed.use = 868)
## PC_ 1 
## Positive:  Apoe, Cst7, Fth1, Fabp5, Ctsb, Cd63, Lyz2, Ctsd, Ctsz, Cd52 
##     Igf1, Lpl, Tyrobp, Ftl1, B2m, Mif, Gnas, Npc2, Cd9, Ank 
##     Csf1, Lgals3bp, Ccl3, Itgax, Hif1a, Eef1a1, Ch25h, Lilrb4a, H2-D1, Ctsl 
## Negative:  Crybb1, Ecscr, Hpgd, Slc2a5, Lrba, Camk2n1, Filip1l, Stab1, Slc40a1, Il7r 
##     Sox4, Lst1, Ddit4, Bank1, Rtn1, Tent5a, Fam102b, Itga9, Serpinf1, Btg2 
##     Fcrls, Lifr, Gp9, Plxna4, Klf7, Upk1b, Cdk5r1, Nav3, Gbp7, Inpp4b 
## PC_ 2 
## Positive:  Ccl12, Iigp1, Ccl2, Irf7, Fcgr4, C4b, Zbp1, Slfn5, Fgl2, Rsad2 
##     Stat1, Slfn2, Gbp2, H2-Q7, Oasl2, Oasl1, H2-K1, Cd72, Bst2, Cxcl10 
##     Rtp4, Tspo, Usp18, Ly6a, Phf11b, H2-D1, Rnf213, Parp14, Ms4a4c, Tor3a 
## Negative:  Trem2, Ctsd, Ank, Tyrobp, Igf1, Ramp1, Serpine2, Pmp22, Ccl6, Hpgd 
##     Mamdc2, Dkk2, Baiap2l2, Hexa, Ang, Cd68, Fabp5, Itgax, Gpnmb, Syngr1 
##     Aplp2, Cadm1, Scd2, Fcrls, Cst7, Cd9, Creg1, Sox4, Arap2, St8sia6 
## PC_ 3 
## Positive:  Stmn1, Pclaf, Spc24, Prc1, Pbk, Tk1, Fcrls, Knl1, Esco2, Lockd 
##     Ccna2, Smc2, Ezh2, Actg1, Ccl4, Spc25, Asf1b, Ccnf, Ran, Racgap1 
##     Calr, Fam102b, Ncl, Lmnb1, Mis18bp1, Ncapg2, Tagln2, Kif15, Slc15a3, Tubb5 
## Negative:  Apoe, Fau, H2-D1, B2m, Lyz2, H2-K1, Fth1, Eef1a1, Eef1b2, H2-Q7 
##     Zfas1, Gas5, Timp2, Uqcrh, Cd52, Rack1, Cd74, Tyrobp, Trem2, Npc2 
##     Uba52, Ctsl, Cd63, Apoc1, C4b, Lgals3bp, Ccl5, H2-Ab1, Ly86, Crlf2 
## PC_ 4 
## Positive:  Serpine2, Grn, Tyrobp, Ramp1, Hexa, Oasl2, Syngr1, Arsb, Trem2, Ctsd 
##     Rtp4, Baiap2l2, Irf7, Gnas, Ly86, Rsad2, Fkbp2, Slfn5, Bst2, Trim30a 
##     Cst7, Slc2a5, Oasl1, Herc6, Ang, Usp18, Psat1, Xaf1, Crybb1, Fam102b 
## Negative:  Lgals3, Lilrb4a, Atf3, Id2, Ccl4, Lilr4b, Fn1, Vim, Apoe, Iqgap1 
##     Mmp12, Bcl2a1d, Cybb, Spp1, Fth1, Glipr1, Plaur, Lpl, Plk2, Emp3 
##     Spp1-EGFP, Rai14, Cd44, C3, Cxcl16, Slamf7, Nes, Tnfrsf12a, Fam20c, Cd74 
## PC_ 5 
## Positive:  Pclaf, Prc1, Spc24, Stmn1, Pbk, Knl1, Fau, Esco2, Lockd, Tk1 
##     Ccna2, Racgap1, Mis18bp1, Ccnf, Fth1, Spc25, Smc2, Ezh2, Trim59, Kif15 
##     Asf1b, Knstrn, Lig1, Dut, Ncapg2, H2afx, Cit, Cenpk, H2-D1, Diaph3 
## Negative:  Rsad2, Slfn5, Usp18, Oasl1, Parp14, Herc6, Oasl2, Irf7, Ccl4, Cxcl10 
##     Slc15a3, Fcrls, Slfn8, Rtp4, Trim30a, Cmpk2, Rnf213, Ccl3, Phf11d, Tnfsf10 
##     Znfx1, Stat2, Ddx58, Atf3, Lilrb4a, C3ar1, Cdk6, Iigp1, Irgm1, Xaf1
length(VariableFeatures(GEX.seur))
## [1] 1390
head(VariableFeatures(GEX.seur),300)
##   [1] "Cxcl10"    "Spp1"      "Ccl4"      "Cxcl13"    "Ccl5"      "Ccl3"     
##   [7] "Cxcl9"     "Ccl12"     "Spp1-EGFP" "Mmp12"     "Cd74"      "H2-Aa"    
##  [13] "Gpnmb"     "Ccl2"      "Rsad2"     "Cst7"      "Cck"       "H2-Ab1"   
##  [19] "Lpl"       "Ttr"       "Lgals3"    "Fabp5"     "Iigp1"     "Fn1"      
##  [25] "H2-Eb1"    "Apoc1"     "Ccl6"      "Gdf15"     "Saa3"      "Apoe"     
##  [31] "Enpp2"     "Wfdc17"    "Ch25h"     "Il1b"      "Prc1"      "Il1rn"    
##  [37] "Lyz2"      "Tnf"       "Ccl7"      "Stmn1"     "Pclaf"     "Cxcl2"    
##  [43] "Egr3"      "Igf1"      "Dqx1"      "Cd5l"      "Aldh1a3"   "Igfbp5"   
##  [49] "Ank"       "Cdkn1a"    "Ly6k"      "Apoc4"     "Slc7a11"   "Wdcp"     
##  [55] "Ms4a7"     "Clec4d"    "Itgax"     "Prdx1"     "Cox6a2"    "Cd52"     
##  [61] "Ccl8"      "Csf1"      "Atf3"      "Ctla2a"    "Ifnb1"     "Ctsd"     
##  [67] "Cd72"      "Syngr1"    "Cybb"      "Gbp2"      "Gnas"      "Vim"      
##  [73] "Ctsb"      "Bst2"      "Ly6a"      "Postn"     "Lilrb4a"   "Tyrobp"   
##  [79] "Mif"       "Cd69"      "Neil3"     "Aplp2"     "Plp1"      "C4b"      
##  [85] "Dennd5b"   "Nr4a1"     "Fth1"      "Edn1"      "Lgals1"    "Ftl1"     
##  [91] "Cd63"      "Olr1"      "Ramp1"     "Pbld1"     "Lgals3bp"  "Chac1"    
##  [97] "Bex1"      "Pbk"       "Slfn5"     "Trib3"     "Clu"       "Knl1"     
## [103] "Hif1a"     "Ccr1"      "Flt1"      "Serpine2"  "Mamdc2"    "Gadd45b"  
## [109] "Ctsl"      "Pou3f1"    "Gria2"     "Sprr1a"    "Oasl1"     "Ccl9"     
## [115] "Colec12"   "Ctsz"      "Tcap"      "Id3"       "Fcgr4"     "Ccdc63"   
## [121] "Ms4a4b"    "Id2"       "Esco2"     "Dkk2"      "Il4i1"     "Gdf3"     
## [127] "Grb14"     "Gpr85"     "Rgs1"      "Irf7"      "Ptgds"     "Fxyd5"    
## [133] "Enah"      "Cstb"      "Lgi2"      "Lockd"     "Phlda3"    "Trem2"    
## [139] "Cxcl14"    "Cd9"       "Acod1"     "Baiap2l2"  "Crip1"     "Ecrg4"    
## [145] "Mt1"       "Oasl2"     "Timp2"     "Ccl22"     "Clec12a"   "Egfl7"    
## [151] "Rhox5"     "Spag5"     "H2-Q7"     "Gapdh"     "Loxl2"     "Slamf7"   
## [157] "Tspo"      "Uchl1"     "Ccna2"     "Stap1"     "Ptprf"     "Ppp1r14a" 
## [163] "Alpk2"     "Crlf2"     "Plac8"     "Aspm"      "Rab3c"     "Ccnb1"    
## [169] "Hbegf"     "Capg"      "Gfod2"     "Usp18"     "Snca"      "Sgo2a"    
## [175] "Mgl2"      "Shcbp1"    "Gpr84"     "Pkm"       "Hp"        "Kcnk9"    
## [181] "Pld3"      "Fabp3"     "Plk2"      "Cspg4"     "Nes"       "Spc24"    
## [187] "Cd40"      "Cd83"      "Kcnq1ot1"  "Scd2"      "Ankle1"    "Slc7a14"  
## [193] "Tlr2"      "Tnip3"     "Rcan1"     "Ms4a4c"    "Tk1"       "Slfn2"    
## [199] "Ptgs2"     "Tbx19"     "Meis2"     "Lrrc31"    "Sh3d19"    "Rgs17"    
## [205] "Olfr889"   "Stat1"     "Peg3"      "Tnfrsf11b" "Ptprz1"    "Axl"      
## [211] "Zbp1"      "Madcam1"   "Lhfpl4"    "Apod"      "H2-D1"     "Adamts1"  
## [217] "Mal"       "Plau"      "Pf4"       "Atp6v0d2"  "Cd68"      "Syt5"     
## [223] "Selenow"   "H2-K1"     "Slc1a2"    "Fam20c"    "Serpina3g" "Cd36"     
## [229] "Cit"       "Iqgap1"    "Smc2"      "Pcgf2"     "Ptgs2os2"  "Smtnl1"   
## [235] "Mir670hg"  "Gchfr"     "Cldn11"    "Sh3gl2"    "Elavl2"    "Osbp2"    
## [241] "Sdk2"      "Fzd3"      "Gpc4"      "Amot"      "Kif28"     "Pdpn"     
## [247] "Pnmal2"    "Gabrb3"    "Olfr56"    "Igsf11"    "Arx"       "Fgl2"     
## [253] "Cxcl16"    "Npc2"      "Naaa"      "Sag"       "Lrpap1"    "Cd38"     
## [259] "Anks1b"    "Serpinb8"  "Ccne1"     "Phf11b"    "Cpd"       "Apoc2"    
## [265] "Bcl2a1d"   "C2"        "Gfap"      "Tpi1"      "Grn"       "Anxa5"    
## [271] "Cp"        "Ptprd"     "Plaur"     "Rtp4"      "Gusb"      "Aspa"     
## [277] "Egr2"      "Nfasc"     "Cd3d"      "Btnl4"     "Kcnq3"     "Crybb1"   
## [283] "Myo1e"     "Lat"       "Ctnna3"    "Syt1"      "H2afz"     "Bhlhe40"  
## [289] "Prr11"     "Parp14"    "Cadm1"     "Ier3"      "Mxd3"      "Hexa"     
## [295] "Hpse"      "St8sia6"   "Ldha"      "Slfn1"     "Gbp4"      "Lilr4b"
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt, ) +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "cnt", cols = color.cnt)

#GEX.seur@reductions$pca@cell.embeddings

check PC_1/2 median/mean

GEX.seur@meta.data[,c("PC_1","PC_2")] <- GEX.seur@reductions$pca@cell.embeddings[rownames(GEX.seur@meta.data),c("PC_1","PC_2")]
cowplot::plot_grid(
as.data.frame(GEX.seur@meta.data) %>% select(cnt,PC_1,PC_2) %>%
  group_by(cnt) %>%
  summarize(m1=mean(PC_1), m2=mean(PC_2),
            n1=median(PC_1), n2=median(PC_2)) %>%
  ggplot(aes(x=m1, y=m2, color=cnt, label=cnt)) + geom_point(color=color.cnt) + geom_text_repel(show.legend = F,color=color.cnt, size=3.5) + labs(title="PC Mean") + theme_classic2(),
as.data.frame(GEX.seur@meta.data) %>% select(cnt,PC_1,PC_2) %>%
  group_by(cnt) %>%
  summarize(m1=mean(PC_1), m2=mean(PC_2),
            n1=median(PC_1), n2=median(PC_2)) %>%
  ggplot(aes(x=n1, y=n2, color=cnt, label=cnt)) + geom_point(color=color.cnt) + geom_text_repel(show.legend = F,color=color.cnt, size=3.5) + labs(title="PC Median") + theme_classic2(),
ncol = 2)

check PCA

cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Spp1", reduction = "pca",dims = 1:2),
  DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt),
rel_widths = c(5,5.6),
ncol = 2)

cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Spp1-EGFP", reduction = "pca",dims = 1:2),
  DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt),
rel_widths = c(5,5.6),
ncol = 2)

cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "P2ry12", reduction = "pca",dims = 1:2),
  DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt),
rel_widths = c(5,5.6),
ncol = 2)

cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Mki67", reduction = "pca",dims = 1:2),
  DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt),
rel_widths = c(5,5.6),
ncol = 2)

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", split.by = "cnt", cols = color.cnt, ncol = 3)

DimHeatmap(GEX.seur, dims = 1:12, cells = 6000, balanced = TRUE,ncol = 4,nfeatures = 80)

(GEX.seur@reductions$pca@feature.loadings %>% as.data.frame)[,1:8] %>% arrange(desc(PC_1) ) %>% head(40)
##                PC_1         PC_2          PC_3          PC_4         PC_5
## Apoe     0.12500859  0.021924692 -0.1441137617 -0.0905009225  0.049730734
## Cst7     0.12238507 -0.048085902 -0.0358791868  0.0679274316  0.007009488
## Fth1     0.11034388  0.042927671 -0.0925807309 -0.0789392314  0.088070464
## Fabp5    0.11016458 -0.052799666 -0.0047849246 -0.0171977019 -0.021476707
## Ctsb     0.10990714 -0.008559184 -0.0281109256  0.0111130622 -0.023351234
## Cd63     0.10818542  0.009186587 -0.0488474404 -0.0351982363  0.027076546
## Lyz2     0.10678830 -0.006881866 -0.1042972513 -0.0183104502  0.054388518
## Ctsd     0.10574859 -0.089161487 -0.0219451625  0.0817485281 -0.030483656
## Ctsz     0.10441218 -0.012486622 -0.0276072630  0.0099890858  0.002926231
## Cd52     0.10227758  0.056533646 -0.0589972148 -0.0063074874  0.065536423
## Igf1     0.10020421 -0.070846684  0.0085575795 -0.0066821365 -0.039675793
## Lpl      0.09915171 -0.041189055  0.0528122116 -0.0752158158 -0.027318480
## Tyrobp   0.09599701 -0.074658856 -0.0572709372  0.0964621518  0.040726996
## Ftl1     0.09472093 -0.014117507 -0.0360902957 -0.0496492969  0.038236325
## B2m      0.09174523  0.078179665 -0.1107401923  0.0286066708  0.066697121
## Mif      0.09034305 -0.012032564  0.0244762234  0.0345747227  0.004733014
## Gnas     0.08854069 -0.029015510  0.0140652656  0.0774564623  0.010335327
## Npc2     0.08581197  0.002578323 -0.0534091808  0.0414969809  0.031483106
## Cd9      0.08562786 -0.047567501 -0.0031681715 -0.0200681472 -0.005125874
## Ank      0.08536142 -0.078449472  0.0191939352  0.0607812243 -0.025127000
## Csf1     0.08531621 -0.030680608  0.0296871408 -0.0236623880 -0.055435412
## Lgals3bp 0.08396246  0.069102467 -0.0468030777  0.0482506858 -0.008417322
## Ccl3     0.08213531 -0.022668987  0.0516302312 -0.0533060478 -0.065083854
## Itgax    0.08206594 -0.052744229 -0.0157853815  0.0025968625 -0.014005518
## Hif1a    0.08068022 -0.044666760  0.0177756171  0.0384092967 -0.045904385
## Eef1a1   0.08055212 -0.006837531 -0.0870921245  0.0077676585  0.054921588
## Ch25h    0.08044329  0.006519790  0.0685956613 -0.0335148396 -0.050724240
## Lilrb4a  0.07924947 -0.009810840  0.0485417523 -0.1272634887 -0.058516403
## H2-D1    0.07916361  0.093661223 -0.1244507098 -0.0456695562  0.069070390
## Ctsl     0.07899560 -0.008789674 -0.0503094611  0.0549571938  0.002892127
## Aplp2    0.07885196 -0.051408165  0.0210250705 -0.0045263253 -0.039743192
## Cox6a2   0.07773751 -0.024552437 -0.0150668438  0.0492204722  0.007773934
## Timp2    0.07754955 -0.021722874 -0.0612060248  0.0313058125  0.013932664
## Spp1     0.07708269 -0.013145039  0.0087505959 -0.0791512053 -0.011659033
## Axl      0.07636212  0.010533011 -0.0245902011  0.0003342218 -0.005580632
## Syngr1   0.07569752 -0.052629108  0.0007126634  0.0835071580  0.005345428
## Gapdh    0.07566385  0.000574007  0.0455580367  0.0347961106  0.015915273
## Pld3     0.07553428 -0.037776831 -0.0096437567 -0.0032568377 -0.017028082
## Fau      0.07479575  0.015691137 -0.1412978100  0.0044090514  0.103813309
## Crlf2    0.07335061 -0.028623250 -0.0398285034  0.0228477374  0.019167524
##                   PC_6          PC_7         PC_8
## Apoe      0.0982195466 -0.0983424473 -0.008502093
## Cst7      0.0290311737 -0.0099327053 -0.002512828
## Fth1     -0.0297440453 -0.0051944329  0.047743260
## Fabp5     0.0310095754 -0.0216505138  0.072161948
## Ctsb     -0.0224101180 -0.0964747600 -0.105253709
## Cd63     -0.0242579948 -0.0271609402 -0.058514049
## Lyz2      0.0540286346 -0.0302003608 -0.032909971
## Ctsd     -0.0381446250 -0.1027168315 -0.083474584
## Ctsz     -0.0297080229 -0.0130895378 -0.091642101
## Cd52     -0.0358904233  0.0140737293 -0.019318791
## Igf1      0.0557677251 -0.0004352220  0.086835198
## Lpl       0.0002024007  0.0422521666  0.060760529
## Tyrobp   -0.0747749449  0.0071451541  0.027927325
## Ftl1     -0.0910036900 -0.0989819743  0.024889730
## B2m      -0.0098065539 -0.0077260605 -0.118524401
## Mif      -0.0799994811  0.0130145225  0.012739368
## Gnas     -0.0612871806 -0.0118103754 -0.029611396
## Npc2     -0.0293213849 -0.0266011510 -0.026496263
## Cd9      -0.0462955169 -0.0418240993 -0.074897457
## Ank       0.0490732929  0.0654246650  0.004755194
## Csf1      0.0448873163  0.0398570863 -0.012255144
## Lgals3bp -0.0096465214 -0.0769078085 -0.103942311
## Ccl3     -0.0278761076 -0.0002478853 -0.055699267
## Itgax     0.0814978539  0.0501108502  0.050126671
## Hif1a     0.0292531506 -0.0112293694 -0.056199937
## Eef1a1   -0.0029653472 -0.0230429748  0.004986999
## Ch25h    -0.0296625436  0.0078588754  0.024366861
## Lilrb4a   0.0068344690 -0.0458841080  0.042066438
## H2-D1     0.0474573253 -0.0031320050 -0.111664276
## Ctsl     -0.0243827936  0.0233000414 -0.086850952
## Aplp2     0.0661102779 -0.0003045849  0.005799693
## Cox6a2    0.0236124750  0.0596410625  0.026435133
## Timp2    -0.0130179956 -0.0894959304 -0.032941522
## Spp1     -0.0041292661 -0.0105701894  0.156801036
## Axl       0.0335174446 -0.0126276303 -0.019725643
## Syngr1   -0.0616502242 -0.0423401940 -0.012947192
## Gapdh    -0.1017024202  0.0255953539  0.003827144
## Pld3      0.0039357945 -0.0613072216 -0.004684004
## Fau      -0.0318015558 -0.0656638937  0.027016024
## Crlf2     0.0437455924  0.0284043632 -0.004665509
(GEX.seur@reductions$pca@feature.loadings %>% as.data.frame)[,1:8] %>% arrange(PC_1 ) %>% head(40)
##                 PC_1          PC_2         PC_3          PC_4         PC_5
## Crybb1   -0.07125275 -0.0174163744  0.061637045  0.0622269759 -0.043559430
## Ecscr    -0.07005891 -0.0191041466  0.064442549  0.0267013693 -0.037679233
## Hpgd     -0.06221909 -0.0597653126  0.062318226  0.0457007202 -0.031627713
## Slc2a5   -0.05589064 -0.0419630764  0.057743705  0.0668063712 -0.040937057
## Lrba     -0.05361008 -0.0325767967  0.057948191  0.0345975740 -0.043381799
## Camk2n1  -0.04435957 -0.0098696682  0.045819352  0.0169471813 -0.026221481
## Filip1l  -0.04377899 -0.0011387692  0.040880296  0.0095577183 -0.035658785
## Stab1    -0.04289914 -0.0222728250  0.056996894  0.0110011530 -0.047525520
## Slc40a1  -0.04203156 -0.0318016407  0.047780477  0.0316343685 -0.047249201
## Il7r     -0.03997281 -0.0196073318  0.041092944  0.0206781712 -0.027762018
## Sox4     -0.03973813 -0.0472777076  0.055569593  0.0256207494 -0.041588757
## Lst1     -0.03859512  0.0260621751 -0.023167824 -0.0087605639  0.019426573
## Ddit4    -0.03662515 -0.0278450020  0.043430091  0.0233451253 -0.033957462
## Bank1    -0.03551595 -0.0221843210  0.049591048  0.0188714562 -0.038604348
## Rtn1     -0.03501590 -0.0174238003  0.032253465  0.0184276814 -0.019850249
## Tent5a   -0.03351556  0.0125236395  0.011491486  0.0212649218 -0.044140608
## Fam102b  -0.03293535 -0.0346810405  0.073080570  0.0619742028 -0.039531649
## Itga9    -0.03226772 -0.0170318533  0.036563965  0.0142177953 -0.034350696
## Serpinf1 -0.03162115 -0.0138954220  0.033631504  0.0078619090 -0.003371705
## Btg2     -0.03103233 -0.0123542271  0.027191907 -0.0147693919 -0.033752866
## Fcrls    -0.03086449 -0.0482267162  0.088309504  0.0591203313 -0.076031969
## Lifr     -0.02856426 -0.0092802571  0.033081122  0.0131790831 -0.028471524
## Gp9      -0.02829243 -0.0176643324  0.033464946  0.0225743866 -0.023541054
## Plxna4   -0.02802836 -0.0132848810  0.032468358  0.0125896755 -0.028106879
## Klf7     -0.02795906 -0.0077095681  0.034270729 -0.0050733047 -0.020162666
## Upk1b    -0.02596451 -0.0283531600  0.046169394  0.0259925490 -0.028292079
## Cdk5r1   -0.02557628 -0.0146448972  0.037416482  0.0238729647 -0.024323067
## Nav3     -0.02448396  0.0223964283  0.040412767  0.0297930532 -0.036329002
## Gbp7     -0.02304218  0.0550766107  0.024410607  0.0230135826 -0.033562917
## Inpp4b   -0.02303367 -0.0210410558  0.030079445  0.0186680503 -0.033392330
## Nfia     -0.02283852 -0.0156165744  0.034132616  0.0205366355 -0.023408068
## Klk8     -0.02155732 -0.0093827933  0.007603048  0.0009380316  0.017973687
## Khdrbs3  -0.02155671 -0.0301695094  0.044366234  0.0350050297 -0.023831469
## Cask     -0.02138948 -0.0174797001  0.031248999  0.0187915466 -0.031402407
## Arid4a   -0.01845172  0.0001199639  0.019577780  0.0071954503 -0.021459080
## Hlf      -0.01843268 -0.0117750565  0.022267571  0.0133176358 -0.014207626
## Fry      -0.01833700 -0.0149468799  0.025258792  0.0253735495 -0.023778229
## Ttc28    -0.01809531 -0.0130445248  0.022546883  0.0124340781 -0.021437365
## Pde4b    -0.01661339 -0.0032311142  0.016503162 -0.0064767670 -0.022405595
## Cdkn1c   -0.01586221 -0.0136016631  0.007934058  0.0071646694 -0.007426902
##                   PC_6          PC_7          PC_8
## Crybb1   -0.0633415123 -3.077838e-02  0.0227968064
## Ecscr    -0.0589302921 -4.044294e-02  0.0095210819
## Hpgd     -0.0225862198  1.823019e-02  0.0466518379
## Slc2a5   -0.0246106177  2.629830e-02  0.0119073416
## Lrba     -0.0278189644 -2.636712e-02  0.0176445691
## Camk2n1  -0.0329895618 -4.933705e-03 -0.0056216252
## Filip1l   0.0107368652  1.123482e-02  0.0095499267
## Stab1    -0.0449203732 -6.037662e-02 -0.0063307953
## Slc40a1   0.0007149646 -4.294400e-02  0.0021210651
## Il7r     -0.0084401180 -1.127658e-02  0.0387445405
## Sox4     -0.0101545132 -2.138098e-02  0.0093627483
## Lst1     -0.0287745399 -2.165201e-02 -0.0225897213
## Ddit4    -0.0176411774 -1.390657e-02  0.0351905738
## Bank1    -0.0105288307 -2.027320e-02  0.0211278451
## Rtn1     -0.0179670528  7.098435e-03  0.0239164333
## Tent5a    0.0204204752 -3.857167e-02  0.0222459453
## Fam102b  -0.0392025699  7.066474e-03 -0.0257295658
## Itga9    -0.0073982446  1.401138e-03  0.0169477878
## Serpinf1 -0.0184724875  6.733560e-03  0.0312529673
## Btg2      0.0287961850 -5.731496e-03 -0.0107138225
## Fcrls    -0.0751583021 -9.902891e-02  0.0048186416
## Lifr      0.0077847703  3.519046e-02  0.0040133212
## Gp9      -0.0248514289  6.246959e-03  0.0294115728
## Plxna4    0.0025864000 -2.539216e-03 -0.0011822520
## Klf7      0.0002998945  1.644622e-02  0.0115026108
## Upk1b    -0.0150634911 -8.489431e-03  0.0370942552
## Cdk5r1   -0.0099464639 -1.941281e-05 -0.0059808086
## Nav3      0.0200091867  8.065609e-02 -0.0596898203
## Gbp7      0.0227850095  5.635619e-02  0.0326802626
## Inpp4b    0.0105538683 -1.359776e-02 -0.0292245987
## Nfia      0.0079608316  7.632060e-03 -0.0083432447
## Klk8     -0.0257461335  4.992017e-02  0.0139343465
## Khdrbs3  -0.0120958931 -7.258875e-03  0.0334096196
## Cask     -0.0120205747 -1.443649e-02 -0.0108235762
## Arid4a    0.0117583299 -1.232380e-02 -0.0054193193
## Hlf      -0.0045868716 -5.792114e-03  0.0193860962
## Fry       0.0197553567 -6.929299e-03  0.0056989181
## Ttc28     0.0128794040 -1.540599e-02 -0.0131753956
## Pde4b    -0.0009610003 -1.690763e-02  0.0009702224
## Cdkn1c   -0.0107060483 -8.448315e-03  0.0159273319
(GEX.seur@reductions$pca@feature.loadings %>% as.data.frame)[,1:8] %>% arrange(desc(PC_2) ) %>% head(40)
##                  PC_1       PC_2          PC_3         PC_4        PC_5
## Ccl12    0.0162794739 0.16742224 -0.0040039504  0.013682355 -0.01507818
## Iigp1    0.0204528910 0.16032868  0.0069663107  0.042628181 -0.05802883
## Ccl2     0.0282896761 0.14263652  0.0520332667  0.002765441 -0.03542026
## Irf7     0.0318651314 0.12331869  0.0049404540  0.077901340 -0.08043553
## Fcgr4    0.0169605197 0.12304679 -0.0195139131 -0.003780755  0.02654228
## C4b      0.0150408778 0.11948059 -0.0480126604 -0.026857729  0.04227467
## Zbp1     0.0347682252 0.11744287  0.0054780670  0.061420669 -0.03491002
## Slfn5    0.0408002900 0.11666886 -0.0113764017  0.073375386 -0.11333964
## Fgl2     0.0266356713 0.11666022 -0.0343098776  0.042119985 -0.01244219
## Rsad2    0.0274427009 0.11661984  0.0155112740  0.074064878 -0.12201823
## Stat1    0.0278988065 0.11350502 -0.0108542966  0.056127469 -0.04931813
## Slfn2    0.0143732799 0.11234312  0.0015865877  0.028124852 -0.04970319
## Gbp2     0.0157322592 0.10708124  0.0144809798  0.020345508 -0.03229534
## H2-Q7    0.0383605171 0.10690828 -0.0677319693 -0.043531962  0.05196830
## Oasl2    0.0399187660 0.10638983 -0.0117524532  0.084516299 -0.08164191
## Oasl1    0.0209166294 0.10600626  0.0121561997  0.064980765 -0.09052601
## H2-K1    0.0683856966 0.10587091 -0.0932349214 -0.032098872  0.04722767
## Cd72     0.0553751986 0.10390699  0.0318429984 -0.062023031  0.02083902
## Bst2     0.0632629709 0.10366641 -0.0177913281  0.069007138 -0.01206454
## Cxcl10   0.0139911607 0.10344971  0.0380972757  0.034365485 -0.07675646
## Rtp4     0.0206657146 0.10173586 -0.0102176914  0.081300333 -0.06916225
## Tspo     0.0542517766 0.09975253 -0.0198969688  0.021404738  0.04125853
## Usp18    0.0267438188 0.09954353 -0.0008885168  0.063381971 -0.09173681
## Ly6a     0.0022563580 0.09798450 -0.0015057321  0.047378693 -0.02565861
## Phf11b   0.0302872845 0.09659267  0.0144178101  0.055511374 -0.05015406
## H2-D1    0.0791636121 0.09366122 -0.1244507098 -0.045669556  0.06907039
## Rnf213   0.0201541702 0.09237321  0.0040935776  0.053486251 -0.06641061
## Parp14   0.0113874444 0.09098651 -0.0020314401  0.060541096 -0.08355942
## Ms4a4c   0.0146097350 0.08601387 -0.0198232618  0.029051633 -0.03635740
## Tor3a    0.0162542058 0.08566886 -0.0023851148  0.023747697 -0.05574633
## Clec2d   0.0058176406 0.08560988 -0.0248405399  0.010669634 -0.03421866
## Trim30a  0.0173357849 0.08496637  0.0031088112  0.068530556 -0.06852049
## Sp100    0.0264298579 0.08314675 -0.0080804610  0.051753760 -0.04999209
## Naaa    -0.0004988087 0.08226191  0.0053919752  0.002839857  0.02936179
## Xaf1     0.0240426349 0.08112075 -0.0036932287  0.063203771 -0.05618673
## Phf11d   0.0201153190 0.08104767  0.0052302414  0.060242881 -0.06134279
## H2-Q6    0.0212337590 0.08095222 -0.0311469889 -0.033554192  0.01989777
## Eif2ak2  0.0186915608 0.07960739  0.0065528684  0.055990260 -0.05175924
## Tap1     0.0172520039 0.07951389  0.0039316491  0.016044144 -0.01639107
## Irgm1    0.0097698549 0.07851849  0.0073718131  0.043003242 -0.05761390
##                 PC_6          PC_7         PC_8
## Ccl12   -0.015831877  0.1064469618 -0.059674320
## Iigp1    0.054274132  0.0523973061  0.099886455
## Ccl2    -0.034584955  0.0754511340  0.004256264
## Irf7    -0.022213058 -0.0623809968  0.041941975
## Fcgr4   -0.008489590  0.1157575343 -0.042743339
## C4b      0.011481987  0.0819919017 -0.072772535
## Zbp1    -0.022081982 -0.0049046570  0.042728270
## Slfn5    0.017816070 -0.0593081461  0.044080833
## Fgl2     0.035680180  0.0590292317 -0.017019890
## Rsad2    0.019488398 -0.1017206572  0.084326537
## Stat1    0.038221939  0.0010588611  0.010306586
## Slfn2    0.015884081 -0.0016120457  0.010574805
## Gbp2     0.018988903  0.0433765500  0.056597300
## H2-Q7    0.027068752  0.0045475046 -0.080500529
## Oasl2    0.020745337 -0.0591439789  0.025343858
## Oasl1    0.007262754 -0.0504890082  0.074773879
## H2-K1    0.023720777  0.0006597087 -0.127019760
## Cd72    -0.076568371  0.0593437795  0.010709725
## Bst2    -0.059694670 -0.0107658271 -0.026332219
## Cxcl10   0.025554260 -0.0393979703  0.067557948
## Rtp4    -0.001607410 -0.0467874242  0.018100120
## Tspo    -0.045782592  0.0707880481  0.026033615
## Usp18    0.033466785 -0.0550425870  0.063094419
## Ly6a    -0.006760161  0.0312121665  0.004691053
## Phf11b  -0.021101275 -0.0220672471  0.047032331
## H2-D1    0.047457325 -0.0031320050 -0.111664276
## Rnf213   0.004368772 -0.0364419249  0.011778933
## Parp14   0.057329511 -0.0386082384  0.007518271
## Ms4a4c   0.016448487 -0.0134554225  0.044999694
## Tor3a    0.014061135 -0.0220619574  0.010343744
## Clec2d   0.041454300 -0.0067255208  0.019840436
## Trim30a  0.022439635 -0.0174939682 -0.005168355
## Sp100    0.016628435 -0.0129979743  0.003116756
## Naaa    -0.006941461  0.1673120028 -0.023255224
## Xaf1     0.009203421 -0.0445425418  0.012291032
## Phf11d   0.011046491 -0.0244776065  0.052817020
## H2-Q6    0.028884681  0.0061324743 -0.066705533
## Eif2ak2  0.023056263 -0.0184201745 -0.004917068
## Tap1     0.023572100 -0.0003693182  0.002300167
## Irgm1    0.020699127 -0.0079206145  0.037388861
(GEX.seur@reductions$pca@feature.loadings %>% as.data.frame)[,1:8] %>% arrange(PC_2) %>% head(40)
##                  PC_1        PC_2          PC_3          PC_4          PC_5
## Trem2     0.072573562 -0.09476902 -0.0535656314  0.0822638219 -0.0006088877
## Ctsd      0.105748592 -0.08916149 -0.0219451625  0.0817485281 -0.0304836564
## Ank       0.085361424 -0.07844947  0.0191939352  0.0607812243 -0.0251270004
## Tyrobp    0.095997007 -0.07465886 -0.0572709372  0.0964621518  0.0407269962
## Igf1      0.100204212 -0.07084668  0.0085575795 -0.0066821365 -0.0396757934
## Ramp1     0.066821945 -0.06859341 -0.0009106848  0.0945544402  0.0039528464
## Serpine2  0.061276155 -0.06241428  0.0177811277  0.1212244665 -0.0423523666
## Pmp22     0.017925858 -0.06227385  0.0530114400  0.0212314553 -0.0432768851
## Ccl6      0.065401480 -0.06169190  0.0239999655  0.0392075675  0.0015770281
## Hpgd     -0.062219085 -0.05976531  0.0623182264  0.0457007202 -0.0316277127
## Mamdc2    0.052140973 -0.05818772 -0.0110836210  0.0370293964 -0.0228383625
## Dkk2      0.050364620 -0.05746608 -0.0044434033  0.0207150685 -0.0291975958
## Baiap2l2  0.068953673 -0.05720752  0.0049583914  0.0810130203 -0.0227623637
## Hexa      0.055474587 -0.05665447 -0.0189513582  0.0891006239 -0.0082415496
## Ang       0.005129050 -0.05642911  0.0132529246  0.0641455992 -0.0120134010
## Cd68      0.064510398 -0.05557507  0.0251678736  0.0203259220 -0.0256276383
## Fabp5     0.110164577 -0.05279967 -0.0047849246 -0.0171977019 -0.0214767070
## Itgax     0.082065945 -0.05274423 -0.0157853815  0.0025968625 -0.0140055181
## Gpnmb     0.063765242 -0.05271795 -0.0055521445 -0.0302818188 -0.0397367580
## Syngr1    0.075697525 -0.05262911  0.0007126634  0.0835071580  0.0053454277
## Aplp2     0.078851963 -0.05140816  0.0210250705 -0.0045263253 -0.0397431921
## Cadm1     0.051739090 -0.05115093  0.0175340291  0.0409213443 -0.0248693470
## Scd2      0.044808079 -0.04901294  0.0372454009  0.0547175707 -0.0246417478
## Fcrls    -0.030864494 -0.04822672  0.0883095035  0.0591203313 -0.0760319691
## Cst7      0.122385071 -0.04808590 -0.0358791868  0.0679274316  0.0070094875
## Cd9       0.085627865 -0.04756750 -0.0031681715 -0.0200681472 -0.0051258738
## Creg1     0.058861859 -0.04751559 -0.0206403713  0.0404787633  0.0107762549
## Sox4     -0.039738127 -0.04727771  0.0555695928  0.0256207494 -0.0415887571
## Arap2     0.037856508 -0.04709568  0.0153970842  0.0351295235 -0.0307207836
## St8sia6   0.044404271 -0.04559247 -0.0095322533  0.0230325811  0.0002929502
## Slc12a2   0.010848179 -0.04544818  0.0200535102  0.0503381735 -0.0407252130
## Crip1     0.064726232 -0.04471265 -0.0111398519  0.0007439775  0.0100403079
## Hif1a     0.080680221 -0.04466676  0.0177756171  0.0384092967 -0.0459043849
## Gpr183   -0.004487438 -0.04395056  0.0650220009 -0.0319637950 -0.0332159172
## Asb10     0.035129675 -0.04362704  0.0049945669  0.0343793328 -0.0272787187
## Slc2a5   -0.055890640 -0.04196308  0.0577437054  0.0668063712 -0.0409370572
## Lgi2      0.042221336 -0.04194324  0.0019697169  0.0607018242 -0.0165784235
## Lpl       0.099151712 -0.04118905  0.0528122116 -0.0752158158 -0.0273184804
## Arsb      0.001852044 -0.04114062  0.0551706173  0.0831666320 -0.0288928962
## Ctsa      0.054584378 -0.03983035  0.0066414219  0.0357924821 -0.0355123625
##                   PC_6          PC_7          PC_8
## Trem2     0.0161269840 -0.0030856035 -3.261265e-02
## Ctsd     -0.0381446250 -0.1027168315 -8.347458e-02
## Ank       0.0490732929  0.0654246650  4.755194e-03
## Tyrobp   -0.0747749449  0.0071451541  2.792733e-02
## Igf1      0.0557677251 -0.0004352220  8.683520e-02
## Ramp1     0.0269497500  0.0269577581  1.337887e-02
## Serpine2  0.0077746157 -0.0234802722 -8.183341e-02
## Pmp22    -0.0037317397 -0.0205957802 -5.736306e-02
## Ccl6     -0.0138745719  0.0710656468  2.586989e-02
## Hpgd     -0.0225862198  0.0182301944  4.665184e-02
## Mamdc2    0.0879639496  0.0490573202  5.848202e-02
## Dkk2      0.0997016520  0.0563821374  8.811760e-02
## Baiap2l2  0.0342400715  0.0304473257 -1.387087e-02
## Hexa     -0.0274126984 -0.0403642148 -4.493049e-02
## Ang      -0.0127478073 -0.0021032115  1.636403e-02
## Cd68     -0.0412620952 -0.0370501503 -2.277190e-02
## Fabp5     0.0310095754 -0.0216505138  7.216195e-02
## Itgax     0.0814978539  0.0501108502  5.012667e-02
## Gpnmb     0.0959186423 -0.0234230321  1.314807e-01
## Syngr1   -0.0616502242 -0.0423401940 -1.294719e-02
## Aplp2     0.0661102779 -0.0003045849  5.799693e-03
## Cadm1     0.0365632596  0.0049759424 -4.812054e-02
## Scd2     -0.0042461769  0.0138118908 -1.427429e-02
## Fcrls    -0.0751583021 -0.0990289054  4.818642e-03
## Cst7      0.0290311737 -0.0099327053 -2.512828e-03
## Cd9      -0.0462955169 -0.0418240993 -7.489746e-02
## Creg1    -0.0149536374 -0.0249466879 -1.488515e-02
## Sox4     -0.0101545132 -0.0213809844  9.362748e-03
## Arap2     0.0508635483  0.0204222622  1.934064e-02
## St8sia6   0.0752622436  0.0647453018  7.530945e-05
## Slc12a2   0.0067656431 -0.0374812828 -8.024335e-02
## Crip1     0.0329516834  0.0189413202  9.111362e-02
## Hif1a     0.0292531506 -0.0112293694 -5.619994e-02
## Gpr183   -0.0023405489  0.0387170331  2.606133e-02
## Asb10     0.0492501521  0.0329847028  9.269871e-03
## Slc2a5   -0.0246106177  0.0262983042  1.190734e-02
## Lgi2      0.0344698563  0.0591095321 -2.163076e-02
## Lpl       0.0002024007  0.0422521666  6.076053e-02
## Arsb     -0.0068382948 -0.0012097709 -1.369025e-02
## Ctsa     -0.0263860297 -0.0812089662 -7.233065e-02
VlnPlot(GEX.seur, features = c("C4b","Fcgr4","Tspo","Naaa","Spp1","Spp1-EGFP"), ncol = 2, group.by = "cnt", pt.size = 0,cols = color.cnt) & 
  geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)& 
    geom_boxplot(outlier.size = 0, fill="white", width=0.15, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=1.2, color="black", alpha=0.55)

decide PCs to use
ElbowPlot(GEX.seur,ndims = 50)

PCs <- 1:18
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 20)
## Computing nearest neighbor graph
## Computing SNN
GEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 1.2)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 21021
## Number of edges: 557343
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7466
## Number of communities: 15
## Elapsed time: 3 seconds

Run UMAP/tSNE

GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 287)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 17:55:01 UMAP embedding parameters a = 0.9922 b = 1.112
## 17:55:01 Read 21021 rows and found 18 numeric columns
## 17:55:01 Using Annoy for neighbor search, n_neighbors = 20
## 17:55:01 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 17:55:04 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\Rtmp08wKfR\file4fe42ef5a37
## 17:55:04 Searching Annoy index using 1 thread, search_k = 2000
## 17:55:09 Annoy recall = 100%
## 17:55:10 Commencing smooth kNN distance calibration using 1 thread
## 17:55:11 Initializing from normalized Laplacian + noise
## 17:55:12 Commencing optimization for 200 epochs, with 617434 positive edges
## 17:55:33 Optimization finished
#saveRDS(GEX.seur,"GEX0829.seur.pure_test.rds")
DimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)

FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

pp.new.1a <- DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
         DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
pp.new.1a

pp.new.1c <- DimPlot(GEX.seur, reduction = "umap", group.by = "seurat_clusters", split.by = "cnt", #cols = color.cnt, 
                     ncol = 3)
pp.new.1c 

pp.new.1d <- DimPlot(GEX.seur, reduction = "umap", group.by = "cnt", split.by = "cnt", cols = color.cnt, ncol = 3)
pp.new.1d 

GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
                                 levels = c(0,1,2,
                                            12,11,
                                            3,7,9,14,
                                            6,10,8,13,
                                            4,5))
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0.05, group.by = "sort_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0, group.by = "sort_clusters")

markers.mig <- c("Ptprc",#"Cd3d","Cd3e","Cd19",
                 "Cd74","Lyz2","Ccl4",
                 "Aif1","P2ry12","C1qa","Spp1",
                 "Top2a","Pcna","Mki67","Mcm6",
                 "Cx3cr1","Il4ra","Il13ra1","Spp1-EGFP",
                 "Fabp5","Hmox1","Ms4a7","Cenpa")
FeaturePlot(GEX.seur, 
            features = markers.mig,
            ncol = 4)

DotPlot(GEX.seur, features = c("Cx3cr1","Spp1","Aif1", # Aif1: Iba1
                               "Fcer1g","Il4ra","Il13ra1"),
        group.by = "sort_clusters")

DotPlot(GEX.seur, features = markers.mig, group.by = "sort_clusters",
                         cols = c("midnightblue","darkorange1")) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)

cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$cnt,
      clusters=GEX.seur$sort_clusters),
                   main = "Cell Count",
                   gaps_row = c(3,6),
      #gaps_col = c(12,13,14),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$cnt,
                                clusters=GEX.seur$sort_clusters)),
                   main = "Cell Ratio",
                   gaps_row = c(3,6),
      #gaps_col = c(12,13,14),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

all markers

# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "sort_clusters"

#GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE, min.pct = 0.05,
#                                  test.use = "MAST",
#                                  #test.use = "wilcox",
#                                  logfc.threshold = 0.25)
GEX.markers.pre <- read.table("sc10x_LYN.marker.pure0829.csv", header = TRUE, sep = ",")
GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
## # A tibble: 120 x 7
## # Groups:   cluster [15]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene   
##        <dbl>      <dbl> <dbl> <dbl>     <dbl>   <int> <chr>  
##  1 0              1.06  0.998 0.91  0               0 P2ry12 
##  2 0              1.01  0.691 0.353 0               0 Crybb1 
##  3 0              0.973 0.819 0.539 0               0 Cd164  
##  4 0              0.895 0.936 0.745 0               0 Maf    
##  5 0              0.880 0.987 0.848 0               0 Fcrls  
##  6 5.28e-309      0.875 0.665 0.347 9.29e-305       0 Ecscr  
##  7 7.70e-249      0.858 0.262 0.061 1.35e-244       0 Ddit4  
##  8 1.02e-241      0.849 0.421 0.172 1.79e-237       0 Slc40a1
##  9 0              0.755 0.993 0.912 0               1 P2ry12 
## 10 8.14e-285      0.607 0.981 0.925 1.43e-280       1 Gpr34  
## # ... with 110 more rows
GEX.markers.pre$cluster <- factor(as.character(GEX.markers.pre$cluster),
                          levels = levels(GEX.seur$sort_clusters))



markers.pre_t32 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.1) %>%
                   top_n(n = 32, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene


markers.pre_t48 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-",GEX.markers.pre$gene,invert = T,value = T)) %>%
                   top_n(n = 48, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t60 <- (GEX.markers.pre %>% group_by(cluster) %>% 
                  filter(pct.1>0.01 & gene %in% grep("Rps|Rpl|mt-",GEX.markers.pre$gene,invert = T,value = T)) %>%
                   top_n(n = 60, wt = avg_log2FC) %>%
                    filter(p_val_adj < 0.01) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene
DotPlot(GEX.seur, features = rev(markers.pre_t60[1:64]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, features = rev(markers.pre_t60[65:128]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, features = rev(markers.pre_t60[129:192]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, features = rev(markers.pre_t60[193:256]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, features = rev(markers.pre_t60[257:320]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, features = rev(markers.pre_t60[321:384]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, features = rev(markers.pre_t60[385:448]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, features = rev(markers.pre_t60[449:508]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

signature score

#### 10x data, calculate signature score       
#
## The code below is from Adam Hamber
## 2D scoring by Itay
get_controls <- function(counts, gene.list, verbose=F, control.genes.per.gene=10)
{
    # Itay: "Such scores are inevitably correlated with cell complexity so to avoid 
    # that I subtract a "control" score which is generated by averaging over a control 
    # gene set. Control gene sets are chosen to contain 100 times more genes than the 
    # real gene set (analogous to averaging over 100 control sets of similar size) and 
    # to have the same distribution of population/bulk - based expression levels as the 
    # real gene set, such that they are expected to have the same number of "zeros" and 
    # to eliminate the correlation with complexity."
    # ---------------------------------------------------------------------------------
    # Going to find control points by finding the closest genes in terms of expression level and % of the time we observe it
    if(verbose){cat(sprintf("Finding %s background genes based on similarity to given gene set [%s genes] \n", 
                            control.genes.per.gene*length(gene.list), length(gene.list)))}
    cat("Summarizing data \n")
    summary = data.frame(gene=row.names(counts), mean.expr = Matrix::rowMeans(counts), fract.zero = Matrix::rowMeans(counts==0), stringsAsFactors = F)
    #summary = data.frame(gene=row.names(counts), mean.expr = apply(counts,1,mean), fract.zero = apply(counts==0,1,mean), stringsAsFactors = F)
    summary$mean.expr.s = scale(summary$mean.expr)
    summary$fract.zero.s = scale(summary$fract.zero)
    actual.genes = summary[summary$gene %in% gene.list,]
    background.genes = summary[!summary$gene %in% gene.list,]
    
    #find the 10 closest genes to each cell cycle marker gene and add them to the lists of control genes
    get_closest_genes <- function(i)
    {
        background.genes$dist = sqrt((background.genes$mean.expr.s - actual.genes$mean.expr.s[i])^2 + 
                                         (background.genes$fract.zero.s - actual.genes$fract.zero.s[i])^2)
        ordered = background.genes$gene[order(background.genes$dist)]
        ordered = ordered[!ordered %in% controls] # don't take genes that already appear in the list 
        closest = head(ordered, n=control.genes.per.gene)
        return(closest)
    }
    controls = c();
    
    for (i in 1:length(gene.list)){
        #info(sprintf("Finding %s control genes for %s", control.genes.per.gene, gene.list[i]))
        closest = get_closest_genes(i)
        #info(sprintf("Found %s: ", length(closest)))
        controls = unique(c(controls, closest))
    }
    
    if(verbose){cat(sprintf("Control gene selection complete. %s genes found. \n", length(controls)))}
    #print(controls)
    return(controls)
}

## Define calculate function
calculate_signature_score <- function(count_matrix, gene_list){
    control_gene <- get_controls(counts = count_matrix,
                                 gene.list = gene_list)
    signature_score <- colMeans(count_matrix[gene_list, ], na.rm = TRUE) - 
        colMeans(count_matrix[control_gene, ], na.rm = TRUE)
    return(signature_score)
}

add_geneset_score <- function(obj, geneset, setname){
  score <- calculate_signature_score(as.data.frame(obj@assays[['RNA']]@data),
                                     geneset)
  obj <- AddMetaData(obj,
                     score,
                     setname)
  return(obj)
}

DAM

DAM.sig <- read.csv("I:/Shared_win/projects/20230811_10x_LYN/analysis_plus_exogene/figures1002/new/ranking_of_DAM_indicator_genes.csv")
colnames(DAM.sig)[1] <- "Ranking"
DAM.sig[1:8,]
##   Ranking     Gene Frequence Total.score
## 1       1     Spp1        11    35.98221
## 2       2     Apoe        11    31.44704
## 3       3   Ifitm3         9    20.82901
## 4       4 Ifi27l2a         9    20.44047
## 5       5     Lyz2        11    20.22463
## 6       6     Cst7         9    19.48372
## 7       7     Cd74         7    18.24720
## 8       8   Lgals3         7    18.24180
DAM.list <- list(top50=DAM.sig$Gene[1:50],
                 top100=DAM.sig$Gene[1:100],
                 top250=DAM.sig$Gene[1:250],
                 top500=DAM.sig$Gene[1:500])
DAM.list
## $top50
##  [1] "Spp1"     "Apoe"     "Ifitm3"   "Ifi27l2a" "Lyz2"     "Cst7"    
##  [7] "Cd74"     "Lgals3"   "Lpl"      "Cd52"     "Ccl5"     "Ccl12"   
## [13] "Lgals3bp" "Cd63"     "H2-Ab1"   "Fn1"      "H2-Aa"    "Fxyd5"   
## [19] "Ccl4"     "Cybb"     "Bst2"     "Cstb"     "H2-Eb1"   "Fth1"    
## [25] "Vim"      "Tspo"     "Ctsb"     "Ccl3"     "Axl"      "Anxa5"   
## [31] "Isg15"    "Lgals1"   "Ccl2"     "Ifi204"   "Igf1"     "Ch25h"   
## [37] "Mif"      "Cxcl10"   "Plin2"    "Fabp5"    "Il1b"     "Crip1"   
## [43] "Slfn2"    "B2m"      "Irf7"     "Cd72"     "Capg"     "Ms4a6c"  
## [49] "Ifit3"    "Apoc1"   
## 
## $top100
##   [1] "Spp1"     "Apoe"     "Ifitm3"   "Ifi27l2a" "Lyz2"     "Cst7"    
##   [7] "Cd74"     "Lgals3"   "Lpl"      "Cd52"     "Ccl5"     "Ccl12"   
##  [13] "Lgals3bp" "Cd63"     "H2-Ab1"   "Fn1"      "H2-Aa"    "Fxyd5"   
##  [19] "Ccl4"     "Cybb"     "Bst2"     "Cstb"     "H2-Eb1"   "Fth1"    
##  [25] "Vim"      "Tspo"     "Ctsb"     "Ccl3"     "Axl"      "Anxa5"   
##  [31] "Isg15"    "Lgals1"   "Ccl2"     "Ifi204"   "Igf1"     "Ch25h"   
##  [37] "Mif"      "Cxcl10"   "Plin2"    "Fabp5"    "Il1b"     "Crip1"   
##  [43] "Slfn2"    "B2m"      "Irf7"     "Cd72"     "Capg"     "Ms4a6c"  
##  [49] "Ifit3"    "Apoc1"    "Fcgr4"    "Il1rn"    "Wfdc17"   "Cxcl2"   
##  [55] "Cxcl16"   "Prdx1"    "Rtp4"     "H2-D1"    "Pkm"      "Stat1"   
##  [61] "Anxa2"    "Gpnmb"    "Zbp1"     "Ftl1"     "Ldha"     "Npc2"    
##  [67] "Ly6a"     "Oasl2"    "Gapdh"    "Prdx5"    "Gbp2"     "Grn"     
##  [73] "Ifi207"   "Ifitm2"   "Tlr2"     "Txn1"     "Phf11b"   "Ctsz"    
##  [79] "Pianp"    "Cd36"     "Itgax"    "Fgl2"     "Ccl6"     "Iigp1"   
##  [85] "Ccl7"     "H2-K1"    "Pld3"     "Tpi1"     "Akr1a1"   "Usp18"   
##  [91] "Rab7b"    "Tmsb10"   "Cd44"     "Aldoa"    "Hcar2"    "Acod1"   
##  [97] "Cd300lb"  "Ctsc"     "Gpr65"    "H2-Q7"   
## 
## $top250
##   [1] "Spp1"     "Apoe"     "Ifitm3"   "Ifi27l2a" "Lyz2"     "Cst7"    
##   [7] "Cd74"     "Lgals3"   "Lpl"      "Cd52"     "Ccl5"     "Ccl12"   
##  [13] "Lgals3bp" "Cd63"     "H2-Ab1"   "Fn1"      "H2-Aa"    "Fxyd5"   
##  [19] "Ccl4"     "Cybb"     "Bst2"     "Cstb"     "H2-Eb1"   "Fth1"    
##  [25] "Vim"      "Tspo"     "Ctsb"     "Ccl3"     "Axl"      "Anxa5"   
##  [31] "Isg15"    "Lgals1"   "Ccl2"     "Ifi204"   "Igf1"     "Ch25h"   
##  [37] "Mif"      "Cxcl10"   "Plin2"    "Fabp5"    "Il1b"     "Crip1"   
##  [43] "Slfn2"    "B2m"      "Irf7"     "Cd72"     "Capg"     "Ms4a6c"  
##  [49] "Ifit3"    "Apoc1"    "Fcgr4"    "Il1rn"    "Wfdc17"   "Cxcl2"   
##  [55] "Cxcl16"   "Prdx1"    "Rtp4"     "H2-D1"    "Pkm"      "Stat1"   
##  [61] "Anxa2"    "Gpnmb"    "Zbp1"     "Ftl1"     "Ldha"     "Npc2"    
##  [67] "Ly6a"     "Oasl2"    "Gapdh"    "Prdx5"    "Gbp2"     "Grn"     
##  [73] "Ifi207"   "Ifitm2"   "Tlr2"     "Txn1"     "Phf11b"   "Ctsz"    
##  [79] "Pianp"    "Cd36"     "Itgax"    "Fgl2"     "Ccl6"     "Iigp1"   
##  [85] "Ccl7"     "H2-K1"    "Pld3"     "Tpi1"     "Akr1a1"   "Usp18"   
##  [91] "Rab7b"    "Tmsb10"   "Cd44"     "Aldoa"    "Hcar2"    "Acod1"   
##  [97] "Cd300lb"  "Ctsc"     "Gpr65"    "H2-Q7"    "Cdkn1a"   "Psat1"   
## [103] "Trim30a"  "Cxcl14"   "Srgn"     "Mmp12"    "Bcl2a1b"  "Tap1"    
## [109] "AB124611" "Xaf1"     "Ly6e"     "Psme1"    "Ctsl"     "Sp100"   
## [115] "Cxcr4"    "Psmb8"    "AA467197" "Pgk1"     "Emp3"     "Csf1"    
## [121] "Mcemp1"   "Gusb"     "Pim1"     "Ctse"     "Cox4i1"   "Il12b"   
## [127] "Msrb1"    "Npm1"     "Alcam"    "Psme2"    "Apoc2"    "Bhlhe40" 
## [133] "Stmn1"    "Gm4951"   "Pla2g7"   "Plaur"    "Tor3a"    "Hspe1"   
## [139] "Tpt1"     "Ndufa1"   "Flt1"     "Tgfbi"    "Cox6c"    "Irgm1"   
## [145] "Ifit2"    "Uba52"    "Igf2r"    "Bola2"    "Ank"      "Tyrobp"  
## [151] "Tpm4"     "Ass1"     "Ms4a4c"   "Ifit1"    "Ybx1"     "Sod2"    
## [157] "Cox8a"    "Fam20c"   "Oas1a"    "Arg1"     "Ms4a7"    "Smpdl3a" 
## [163] "Sh3pxd2b" "Fau"      "Gnas"     "Phf11d"   "Ehd1"     "Saa3"    
## [169] "Cox5a"    "Atox1"    "Id3"      "Lrpap1"   "Olr1"     "Sh3bgrl3"
## [175] "Sash1"    "Hint1"    "Il2rg"    "Ctsd"     "Postn"    "H2-T23"  
## [181] "Ucp2"     "Sdc3"     "Uqcrq"    "Cox6a1"   "Cd300lf"  "Syngr1"  
## [187] "Timp2"    "Atp5e"    "Id2"      "S100a6"   "Cd14"     "Tubb6"   
## [193] "Anp32b"   "Fcgr2b"   "Cd83"     "Psmb9"    "Bcl2a1a"  "Aprt"    
## [199] "Mfsd12"   "Psap"     "Cox6b1"   "Lilr4b"   "Plac8"    "Glrx"    
## [205] "Scimp"    "Rilpl2"   "Psmb6"    "Gm11808"  "Chmp4b"   "Actr3b"  
## [211] "Ly86"     "Fundc2"   "Ifi211"   "Hif1a"    "Serf2"    "Dram1"   
## [217] "Parp14"   "Ptgs2"    "Cxcl9"    "Myo5a"    "Gabarap"  "Cyp4f18" 
## [223] "Shisa5"   "Lilrb4a"  "Cpd"      "Iqgap1"   "Slfn5"    "Tnfaip2" 
## [229] "Nme1"     "Cotl1"    "Tagln2"   "Mt1"      "Mpeg1"    "C3"      
## [235] "Ube2l6"   "Clic4"    "Naaa"     "Gas7"     "Sdcbp"    "Nampt"   
## [241] "Ell2"     "Samhd1"   "Rtcb"     "Eef1g"    "Hmgb2"    "Gng5"    
## [247] "Nfil3"    "Adora1"   "Tpd52"    "Ptger4"  
## 
## $top500
##   [1] "Spp1"          "Apoe"          "Ifitm3"        "Ifi27l2a"     
##   [5] "Lyz2"          "Cst7"          "Cd74"          "Lgals3"       
##   [9] "Lpl"           "Cd52"          "Ccl5"          "Ccl12"        
##  [13] "Lgals3bp"      "Cd63"          "H2-Ab1"        "Fn1"          
##  [17] "H2-Aa"         "Fxyd5"         "Ccl4"          "Cybb"         
##  [21] "Bst2"          "Cstb"          "H2-Eb1"        "Fth1"         
##  [25] "Vim"           "Tspo"          "Ctsb"          "Ccl3"         
##  [29] "Axl"           "Anxa5"         "Isg15"         "Lgals1"       
##  [33] "Ccl2"          "Ifi204"        "Igf1"          "Ch25h"        
##  [37] "Mif"           "Cxcl10"        "Plin2"         "Fabp5"        
##  [41] "Il1b"          "Crip1"         "Slfn2"         "B2m"          
##  [45] "Irf7"          "Cd72"          "Capg"          "Ms4a6c"       
##  [49] "Ifit3"         "Apoc1"         "Fcgr4"         "Il1rn"        
##  [53] "Wfdc17"        "Cxcl2"         "Cxcl16"        "Prdx1"        
##  [57] "Rtp4"          "H2-D1"         "Pkm"           "Stat1"        
##  [61] "Anxa2"         "Gpnmb"         "Zbp1"          "Ftl1"         
##  [65] "Ldha"          "Npc2"          "Ly6a"          "Oasl2"        
##  [69] "Gapdh"         "Prdx5"         "Gbp2"          "Grn"          
##  [73] "Ifi207"        "Ifitm2"        "Tlr2"          "Txn1"         
##  [77] "Phf11b"        "Ctsz"          "Pianp"         "Cd36"         
##  [81] "Itgax"         "Fgl2"          "Ccl6"          "Iigp1"        
##  [85] "Ccl7"          "H2-K1"         "Pld3"          "Tpi1"         
##  [89] "Akr1a1"        "Usp18"         "Rab7b"         "Tmsb10"       
##  [93] "Cd44"          "Aldoa"         "Hcar2"         "Acod1"        
##  [97] "Cd300lb"       "Ctsc"          "Gpr65"         "H2-Q7"        
## [101] "Cdkn1a"        "Psat1"         "Trim30a"       "Cxcl14"       
## [105] "Srgn"          "Mmp12"         "Bcl2a1b"       "Tap1"         
## [109] "AB124611"      "Xaf1"          "Ly6e"          "Psme1"        
## [113] "Ctsl"          "Sp100"         "Cxcr4"         "Psmb8"        
## [117] "AA467197"      "Pgk1"          "Emp3"          "Csf1"         
## [121] "Mcemp1"        "Gusb"          "Pim1"          "Ctse"         
## [125] "Cox4i1"        "Il12b"         "Msrb1"         "Npm1"         
## [129] "Alcam"         "Psme2"         "Apoc2"         "Bhlhe40"      
## [133] "Stmn1"         "Gm4951"        "Pla2g7"        "Plaur"        
## [137] "Tor3a"         "Hspe1"         "Tpt1"          "Ndufa1"       
## [141] "Flt1"          "Tgfbi"         "Cox6c"         "Irgm1"        
## [145] "Ifit2"         "Uba52"         "Igf2r"         "Bola2"        
## [149] "Ank"           "Tyrobp"        "Tpm4"          "Ass1"         
## [153] "Ms4a4c"        "Ifit1"         "Ybx1"          "Sod2"         
## [157] "Cox8a"         "Fam20c"        "Oas1a"         "Arg1"         
## [161] "Ms4a7"         "Smpdl3a"       "Sh3pxd2b"      "Fau"          
## [165] "Gnas"          "Phf11d"        "Ehd1"          "Saa3"         
## [169] "Cox5a"         "Atox1"         "Id3"           "Lrpap1"       
## [173] "Olr1"          "Sh3bgrl3"      "Sash1"         "Hint1"        
## [177] "Il2rg"         "Ctsd"          "Postn"         "H2-T23"       
## [181] "Ucp2"          "Sdc3"          "Uqcrq"         "Cox6a1"       
## [185] "Cd300lf"       "Syngr1"        "Timp2"         "Atp5e"        
## [189] "Id2"           "S100a6"        "Cd14"          "Tubb6"        
## [193] "Anp32b"        "Fcgr2b"        "Cd83"          "Psmb9"        
## [197] "Bcl2a1a"       "Aprt"          "Mfsd12"        "Psap"         
## [201] "Cox6b1"        "Lilr4b"        "Plac8"         "Glrx"         
## [205] "Scimp"         "Rilpl2"        "Psmb6"         "Gm11808"      
## [209] "Chmp4b"        "Actr3b"        "Ly86"          "Fundc2"       
## [213] "Ifi211"        "Hif1a"         "Serf2"         "Dram1"        
## [217] "Parp14"        "Ptgs2"         "Cxcl9"         "Myo5a"        
## [221] "Gabarap"       "Cyp4f18"       "Shisa5"        "Lilrb4a"      
## [225] "Cpd"           "Iqgap1"        "Slfn5"         "Tnfaip2"      
## [229] "Nme1"          "Cotl1"         "Tagln2"        "Mt1"          
## [233] "Mpeg1"         "C3"            "Ube2l6"        "Clic4"        
## [237] "Naaa"          "Gas7"          "Sdcbp"         "Nampt"        
## [241] "Ell2"          "Samhd1"        "Rtcb"          "Eef1g"        
## [245] "Hmgb2"         "Gng5"          "Nfil3"         "Adora1"       
## [249] "Tpd52"         "Ptger4"        "Eif2ak2"       "Areg"         
## [253] "Rsad2"         "Slc31a1"       "Gm2000"        "Cox7c"        
## [257] "Tmem256"       "Cox5b"         "Cyba"          "Il18bp"       
## [261] "Selenow"       "Myl6"          "H2-DMa"        "Rai14"        
## [265] "Dab2"          "Hmox1"         "Tmsb4x"        "Ndufa2"       
## [269] "Cd68"          "Ccdc86"        "Atp6v1a"       "Cox7a2"       
## [273] "Calm1"         "Uqcrh"         "Socs2"         "Gpi1"         
## [277] "Ubl5"          "Colec12"       "Ifit3b"        "Scpep1"       
## [281] "S100a11"       "F3"            "Ms4a6d"        "Hacd2"        
## [285] "Chil3"         "Tuba1c"        "Rap2b"         "Gp49a"        
## [289] "Milr1"         "Fcgr1"         "Stat2"         "Atp5j2"       
## [293] "Uqcr10"        "Ssr4"          "Bcar3"         "Gm49368"      
## [297] "Tmem106a"      "Tubb5"         "Naca"          "Hspa8"        
## [301] "Atp5k"         "Bag1"          "Sec61b"        "Gyg"          
## [305] "Cox7b"         "Ly6c2"         "Top2a"         "Aldh2"        
## [309] "Dpp7"          "Eif3k"         "Cd48"          "C4b"          
## [313] "Pycard"        "Atp5l"         "Uqcr11"        "Osm"          
## [317] "Prelid1"       "Rnf213"        "Prdx4"         "Arpc1b"       
## [321] "Ndufv3"        "Sp110"         "Ndufa13"       "Abca1"        
## [325] "Gm1673"        "Wdr89"         "Sh3bgrl"       "Smdt1"        
## [329] "Gatm"          "Gpr84"         "Slamf8"        "Ccrl2"        
## [333] "Tomm7"         "Gas8"          "Ly6i"          "Bcl2a1d"      
## [337] "Esd"           "Dhx58"         "Ctsa"          "Rxrg"         
## [341] "Prdx6"         "Ndufc1"        "Polr2f"        "Sdc4"         
## [345] "Atp5j"         "Litaf"         "Atp6v0e"       "Cspg4"        
## [349] "Ranbp1"        "Ifi35"         "Fcer1g"        "Calm3"        
## [353] "Rhoc"          "Pde4b"         "Atp5g1"        "Gpx3"         
## [357] "Psmb1"         "Gpx1"          "Eef1a1"        "Ly9"          
## [361] "Igtp"          "H2-Q6"         "Herc6"         "Cd9"          
## [365] "Ran"           "Hebp1"         "Eno1"          "Cdk18"        
## [369] "Eef1b2"        "Gbp7"          "Glipr1"        "Atp6v1f"      
## [373] "H2-DMb1"       "Btf3"          "Slc25a3"       "Brdt"         
## [377] "Csf2ra"        "Eif3f"         "Hpse"          "Gm14305"      
## [381] "Gm14410"       "H2-Q4"         "Ndufb9"        "Epsti1"       
## [385] "Dnaaf3"        "Pf4"           "S100a4"        "Atp5h"        
## [389] "Apobec3"       "Hsp90ab1"      "Tnf"           "Ctss"         
## [393] "Gas6"          "Gbp3"          "Gng12"         "Nceh1"        
## [397] "Mgst1"         "Cfl1"          "Dtnbp1"        "Slc2a6"       
## [401] "Eif5a"         "Atp5c1"        "ptchd1"        "Ptchd1"       
## [405] "Psma7"         "Lap3"          "Rbm3"          "Cycs"         
## [409] "Cox6a2"        "Abcg1"         "Prr15"         "Ahnak"        
## [413] "Ndufb7"        "Nrp1"          "Lmnb1"         "Ninj1"        
## [417] "Mrpl33"        "Arpc3"         "Phlda1"        "Acp5"         
## [421] "Atp5g3"        "C5ar1"         "Arl5c"         "Parp9"        
## [425] "Ndufb8"        "Txndc17"       "Tbca"          "Gm49339"      
## [429] "Tma7"          "Fblim1"        "Eif3h"         "Psme2b"       
## [433] "Mrpl30"        "Arpc2"         "Ptma"          "Rac2"         
## [437] "Ctsh"          "Dcstamp"       "Clec4n"        "Dbi"          
## [441] "H2-T22"        "Sem1"          "Msr1"          "Bax"          
## [445] "S100a10"       "Fabp3"         "Snrpg"         "Bcl2a1"       
## [449] "Serpine2"      "Snrpd2"        "Cd180"         "Pgam1"        
## [453] "2310061I04Rik" "S100a1"        "Serpinb6a"     "Actg1"        
## [457] "Uqcc2"         "Tuba1b"        "Neurl2"        "Gcnt2"        
## [461] "Smc2"          "Atp5b"         "Ifih1"         "Nap1l1"       
## [465] "Cd274"         "Rgs1"          "Parp12"        "Serpine1"     
## [469] "Nsa2"          "Cebpb"         "Csf2rb"        "Cfb"          
## [473] "Ndufa6"        "Sdf2l1"        "Anxa4"         "Psma5"        
## [477] "Nfkbiz"        "Ctla2a"        "Atp5o"         "Ube2l3"       
## [481] "Lipa"          "Pfn1"          "Ndufa7"        "Ndufs6"       
## [485] "Psmb10"        "Mapkapk2"      "Aif1"          "Smc4"         
## [489] "Itgb1bp1"      "Nme2"          "Pfdn5"         "Rassf3"       
## [493] "1810058I24Rik" "Pgls"          "Clec4d"        "Mrpl54"       
## [497] "Tmem160"       "Pomp"          "Slc7a11"       "F10"
lapply(DAM.list, length)
## $top50
## [1] 50
## 
## $top100
## [1] 100
## 
## $top250
## [1] 250
## 
## $top500
## [1] 500
GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top50, "DAM.sig_top50")
## Summarizing data
GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top100, "DAM.sig_top100")
## Summarizing data
GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top250, "DAM.sig_top250")
## Summarizing data
GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top500, "DAM.sig_top500")
## Summarizing data
ppnew.2b <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
                        raster = T, pt.size = 3.5,
                        cols = c("lightgrey","red"))#,
            #keep.scale = "all")
ppnew.2b

mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
ppnew.2d <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
                        raster = T, pt.size = 3.5,
            keep.scale = "all") & scale_color_gradientn(colors = rev(mapal))
ppnew.2d

ppnew.2.v1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.2.v1

ppnew.2.v2 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      ncol = 2, pt.size = 0, raster = F) & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=1.6, color="black", alpha=0.55)
ppnew.2.v2

ppnew.3.v1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.3.v1

ppnew.3.v2 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      same.y.lims = F,
                      group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P05.CTR","P05.GFP"),
                                                  c("P28.CTR","P28.GFP")),
                               label.y = c(1,1.4),
                               size=3
                               )
ppnew.3.v2

ppnew.3.v3 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      same.y.lims = F,
                      group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & ylim(c(-0.55,1.55)) &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("EAE.CTL","EAE.MIG"),
                                                  c("EAE.MIG","EAE.GFP"),
                                                  c("EAE.CTL","EAE.GFP"),
                                                  c("SIM.CTL","SIM.MIG"),
                                                  c("SIM.MIG","SIM.GFP"),
                                                  c("SIM.CTL","SIM.GFP"),
                                                  c("AD.CTL","AD.MIG"),
                                                  c("AD.MIG","AD.GFP"),
                                                  c("AD.CTL","AD.GFP")),
                               label.y = c(0.8,1.15,1.5,0.8,1.15,1.5,0.8,1.15,1.5),
                               size=1.5
                               )
ppnew.3.v3

module

module.path <- "I:/Shared_win/projects/20230211_microglia_module/analysis_new202310/result.new/modules.GO/"

module.list <- list(m1.blue=as.vector(unlist(read.table(paste0(module.path,"module01.blue.txt")))),
                 m2.purple=as.vector(unlist(read.table(paste0(module.path,"module02.purple.txt")))),
                 m3.darkorange=as.vector(unlist(read.table(paste0(module.path,"module03.darkorange.txt")))),
                 m4.yellow=as.vector(unlist(read.table(paste0(module.path,"module04.yellow.txt")))),
                 m5.yellowgreen=as.vector(unlist(read.table(paste0(module.path,"module05.yellowgreen.txt")))),
                 m6.green=as.vector(unlist(read.table(paste0(module.path,"module06.green.txt")))),
                 m7.darkgreen=as.vector(unlist(read.table(paste0(module.path,"module07.darkgreen.txt")))),
                 m8.turquoise=as.vector(unlist(read.table(paste0(module.path,"module08.turquoise.txt"))))
                 )
lapply(module.list, length)
## $m1.blue
## [1] 3070
## 
## $m2.purple
## [1] 1286
## 
## $m3.darkorange
## [1] 108
## 
## $m4.yellow
## [1] 69
## 
## $m5.yellowgreen
## [1] 45
## 
## $m6.green
## [1] 375
## 
## $m7.darkgreen
## [1] 700
## 
## $m8.turquoise
## [1] 487
names(module.list)
## [1] "m1.blue"        "m2.purple"      "m3.darkorange"  "m4.yellow"     
## [5] "m5.yellowgreen" "m6.green"       "m7.darkgreen"   "m8.turquoise"
for(nn in names(module.list)){
  GEX.seur <- add_geneset_score(GEX.seur, module.list[[nn]], nn)
}
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data 
## Summarizing data
ppnew.5b <- FeaturePlot(GEX.seur,features = c(names(module.list)), ncol = 4,
                        raster = T, pt.size = 3.5,
                        cols = c("lightgrey","red"))#,
            #keep.scale = "all")
ppnew.5b

ppnew.5d <- FeaturePlot(GEX.seur,features = c(names(module.list)), ncol = 4,
                        raster = T, pt.size = 3.5,
            keep.scale = "all") & scale_color_gradientn(colors = rev(mapal))
ppnew.5d

ppnew.6.v1 <- VlnPlot(GEX.seur, features = c(names(module.list)), 
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.6.v1

ppnew.6.v2 <- VlnPlot(GEX.seur, features = c(names(module.list)), 
                      ncol = 2, pt.size = 0, raster = F) & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=1.6, color="black", alpha=0.55)
ppnew.6.v2

ppnew.7.v1 <- VlnPlot(GEX.seur, features = c(names(module.list)), 
                      group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.7.v1

ppnew.7.v2 <- VlnPlot(GEX.seur, features = c(names(module.list)), 
                      same.y.lims = F,
                      group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P05.CTR","P05.GFP"),
                                                  c("P28.CTR","P28.GFP")),
                               label.y = c(1,1.4),
                               size=3
                               )
ppnew.7.v2

ppnew.7.v3 <- VlnPlot(GEX.seur, features = c(names(module.list)), 
                      same.y.lims = F,
                      group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & ylim(c(-0.2,1.25)) &
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("EAE.CTL","EAE.MIG"),
                                                  c("EAE.MIG","EAE.GFP"),
                                                  c("EAE.CTL","EAE.GFP"),
                                                  c("SIM.CTL","SIM.MIG"),
                                                  c("SIM.MIG","SIM.GFP"),
                                                  c("SIM.CTL","SIM.GFP"),
                                                  c("AD.CTL","AD.MIG"),
                                                  c("AD.MIG","AD.GFP"),
                                                  c("AD.CTL","AD.GFP"),
                                                  c("EAE.GFP","SIM.GFP"),
                                                  c("SIM.GFP","AD.GFP"),
                                                  c("EAE.GFP","AD.GFP")),
                               label.y = c(0.5,0.6,0.7,0.5,0.6,0.7,0.5,0.6,0.7,0.8,0.9,1.0),
                               size=2
                               )
ppnew.7.v3

final features

final.list <- list(a=c("H2-Aa","Cd74","H2-Ab1",
                       "H2-Eb1","Iigp1","Gpx3"),
                   b=c("Spp1","Fabp5","Tspo",
                       "Axl","Ly9"),
                   c=c("Ccl5","Ccl7","Ccl12",
                       "Cxcl10","Stat1"),
                   d=c("Trem2","Grn","Cd9",
                       "Tyrobp","Atp6ap1","Fam20c"),
                   e=c("Tmem119","Cx3cr1","Slc2a5",
                       "Sall1","Tgfbr1"),
                   f=c("Igf1","Gnas","Fdps",
                       "Dhcr7","Cyp51"))
lapply(final.list, function(x){x %in% rownames(GEX.seur)})
## $a
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
## 
## $b
## [1] TRUE TRUE TRUE TRUE TRUE
## 
## $c
## [1] TRUE TRUE TRUE TRUE TRUE
## 
## $d
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
## 
## $e
## [1] TRUE TRUE TRUE TRUE TRUE
## 
## $f
## [1] TRUE TRUE TRUE TRUE TRUE
lapply(final.list, function(x){
  lapply(module.list, function(y){
    x %in% y
  })
  })
## $a
## $a$m1.blue
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $a$m2.purple
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $a$m3.darkorange
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $a$m4.yellow
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $a$m5.yellowgreen
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
## 
## $a$m6.green
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $a$m7.darkgreen
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $a$m8.turquoise
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## 
## $b
## $b$m1.blue
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $b$m2.purple
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $b$m3.darkorange
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $b$m4.yellow
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $b$m5.yellowgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $b$m6.green
## [1] TRUE TRUE TRUE TRUE TRUE
## 
## $b$m7.darkgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $b$m8.turquoise
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## 
## $c
## $c$m1.blue
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $c$m2.purple
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $c$m3.darkorange
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $c$m4.yellow
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $c$m5.yellowgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $c$m6.green
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $c$m7.darkgreen
## [1] TRUE TRUE TRUE TRUE TRUE
## 
## $c$m8.turquoise
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## 
## $d
## $d$m1.blue
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $d$m2.purple
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $d$m3.darkorange
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
## 
## $d$m4.yellow
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $d$m5.yellowgreen
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $d$m6.green
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $d$m7.darkgreen
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $d$m8.turquoise
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## 
## $e
## $e$m1.blue
## [1] FALSE FALSE FALSE FALSE  TRUE
## 
## $e$m2.purple
## [1]  TRUE  TRUE  TRUE  TRUE FALSE
## 
## $e$m3.darkorange
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $e$m4.yellow
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $e$m5.yellowgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $e$m6.green
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $e$m7.darkgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $e$m8.turquoise
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## 
## $f
## $f$m1.blue
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $f$m2.purple
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $f$m3.darkorange
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $f$m4.yellow
## [1] TRUE TRUE TRUE TRUE TRUE
## 
## $f$m5.yellowgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $f$m6.green
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $f$m7.darkgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $f$m8.turquoise
## [1] FALSE FALSE FALSE FALSE FALSE
final.list.rename <- list(m5.yellowgreen=c("H2-Aa","Cd74","H2-Ab1",
                       "H2-Eb1","Iigp1","Gpx3"),
                   m6.green=c("Spp1","Fabp5","Tspo",
                       "Axl","Ly9"),
                   m7.darkgreen=c("Ccl5","Ccl7","Ccl12",
                       "Cxcl10","Stat1"),
                   m3.darkorange=c("Trem2","Grn","Cd9",
                       "Tyrobp","Atp6ap1","Fam20c"),
                   m21.purpleblue=c("Tmem119","Cx3cr1","Slc2a5",
                       "Sall1","Tgfbr1"),
                   m4.yellow=c("Igf1","Gnas","Fdps",
                       "Dhcr7","Cyp51"))
lapply(final.list.rename, function(x){
  lapply(module.list, function(y){
    x %in% y
  })
  })
## $m5.yellowgreen
## $m5.yellowgreen$m1.blue
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $m5.yellowgreen$m2.purple
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $m5.yellowgreen$m3.darkorange
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $m5.yellowgreen$m4.yellow
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $m5.yellowgreen$m5.yellowgreen
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
## 
## $m5.yellowgreen$m6.green
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $m5.yellowgreen$m7.darkgreen
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $m5.yellowgreen$m8.turquoise
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## 
## $m6.green
## $m6.green$m1.blue
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m6.green$m2.purple
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m6.green$m3.darkorange
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m6.green$m4.yellow
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m6.green$m5.yellowgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m6.green$m6.green
## [1] TRUE TRUE TRUE TRUE TRUE
## 
## $m6.green$m7.darkgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m6.green$m8.turquoise
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## 
## $m7.darkgreen
## $m7.darkgreen$m1.blue
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m7.darkgreen$m2.purple
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m7.darkgreen$m3.darkorange
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m7.darkgreen$m4.yellow
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m7.darkgreen$m5.yellowgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m7.darkgreen$m6.green
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m7.darkgreen$m7.darkgreen
## [1] TRUE TRUE TRUE TRUE TRUE
## 
## $m7.darkgreen$m8.turquoise
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## 
## $m3.darkorange
## $m3.darkorange$m1.blue
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $m3.darkorange$m2.purple
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $m3.darkorange$m3.darkorange
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
## 
## $m3.darkorange$m4.yellow
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $m3.darkorange$m5.yellowgreen
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $m3.darkorange$m6.green
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $m3.darkorange$m7.darkgreen
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## $m3.darkorange$m8.turquoise
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
## 
## 
## $m21.purpleblue
## $m21.purpleblue$m1.blue
## [1] FALSE FALSE FALSE FALSE  TRUE
## 
## $m21.purpleblue$m2.purple
## [1]  TRUE  TRUE  TRUE  TRUE FALSE
## 
## $m21.purpleblue$m3.darkorange
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m21.purpleblue$m4.yellow
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m21.purpleblue$m5.yellowgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m21.purpleblue$m6.green
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m21.purpleblue$m7.darkgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m21.purpleblue$m8.turquoise
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## 
## $m4.yellow
## $m4.yellow$m1.blue
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m4.yellow$m2.purple
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m4.yellow$m3.darkorange
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m4.yellow$m4.yellow
## [1] TRUE TRUE TRUE TRUE TRUE
## 
## $m4.yellow$m5.yellowgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m4.yellow$m6.green
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m4.yellow$m7.darkgreen
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $m4.yellow$m8.turquoise
## [1] FALSE FALSE FALSE FALSE FALSE
final.list.rename
## $m5.yellowgreen
## [1] "H2-Aa"  "Cd74"   "H2-Ab1" "H2-Eb1" "Iigp1"  "Gpx3"  
## 
## $m6.green
## [1] "Spp1"  "Fabp5" "Tspo"  "Axl"   "Ly9"  
## 
## $m7.darkgreen
## [1] "Ccl5"   "Ccl7"   "Ccl12"  "Cxcl10" "Stat1" 
## 
## $m3.darkorange
## [1] "Trem2"   "Grn"     "Cd9"     "Tyrobp"  "Atp6ap1" "Fam20c" 
## 
## $m21.purpleblue
## [1] "Tmem119" "Cx3cr1"  "Slc2a5"  "Sall1"   "Tgfbr1" 
## 
## $m4.yellow
## [1] "Igf1"  "Gnas"  "Fdps"  "Dhcr7" "Cyp51"
lapply(final.list.rename,function(x){
  FeaturePlot(GEX.seur,features = c(x), ncol = 3,
                        raster = T, pt.size = 3.5,
                        cols = c("lightgrey","red"))
})
## $m5.yellowgreen

## 
## $m6.green

## 
## $m7.darkgreen

## 
## $m3.darkorange

## 
## $m21.purpleblue

## 
## $m4.yellow

lapply(final.list.rename,function(x){
  FeaturePlot(GEX.seur,features = c(x), ncol = 3,
                        raster = T, pt.size = 3.5,
                        cols = c("lightgrey","red"))& scale_color_gradientn(colors = rev(mapal))
})
## $m5.yellowgreen

## 
## $m6.green

## 
## $m7.darkgreen

## 
## $m3.darkorange

## 
## $m21.purpleblue

## 
## $m4.yellow

#saveRDS(GEX.seur,"./GEX0829.seur.pure_sort.rds")

forGEO

GEX.seur <- readRDS("./GEX0829.seur.pure_sort.rds")
GEX.seur
## An object of class Seurat 
## 17590 features across 21021 samples within 1 assay 
## Active assay: RNA (17590 features, 1390 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
GEX.seur@meta.data[,grep("snn|pANN",colnames(GEX.seur@meta.data),value = T)] <- NULL
head(GEX.seur@meta.data)
##                    orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAGTCCG-1     Spp1ND       1777          950   2.138436   8.610017
## AAACCCAAGAGAAGGT-1     Spp1ND       2811         1304   1.529705  12.877979
## AAACCCAAGCACTCAT-1     Spp1ND       3211         1426   3.083152  15.789474
## AAACCCAAGTCATCCA-1     Spp1ND       3553         1643   2.420490  11.004785
## AAACCCAAGTCGGCCT-1     Spp1ND       6412         1938   1.325639  17.358079
## AAACCCACAAGAGGTC-1     Spp1ND       3014         1054   2.687459  20.902455
##                     FB.info      S.Score   G2M.Score Phase seurat_clusters
## AAACCCAAGAAGTCCG-1  SIM.MIG -0.001773983 -0.00433029    G1               1
## AAACCCAAGAGAAGGT-1 AD.GFP.1 -0.034564958  0.01326095   G2M               4
## AAACCCAAGCACTCAT-1 AD.GFP.1  0.004795299 -0.01011902     S               9
## AAACCCAAGTCATCCA-1  SIM.CTL  0.003007456 -0.01007400     S               1
## AAACCCAAGTCGGCCT-1  EAE.GFP  0.076364443 -0.06412611     S               9
## AAACCCACAAGAGGTC-1  EAE.GFP -0.010117249 -0.01113632    G1               7
##                        cnt sort_clusters DoubletFinder0.05 DoubletFinder0.1
## AAACCCAAGAAGTCCG-1 SIM.MIG             1           Singlet          Singlet
## AAACCCAAGAGAAGGT-1  AD.GFP             4           Singlet          Singlet
## AAACCCAAGCACTCAT-1  AD.GFP             9           Singlet          Singlet
## AAACCCAAGTCATCCA-1 SIM.CTL             1           Singlet          Singlet
## AAACCCAAGTCGGCCT-1 EAE.GFP             9           Singlet          Singlet
## AAACCCACAAGAGGTC-1 EAE.GFP             7           Singlet          Singlet
##                      preAnno       PC_1        PC_2 DAM.sig_top50
## AAACCCAAGAAGTCCG-1 Microglia -5.4471120 -1.41232903   -0.12886958
## AAACCCAAGAGAAGGT-1 Microglia  0.7139283 -2.41197770    0.05988822
## AAACCCAAGCACTCAT-1 Microglia  6.7730504 -1.69563381    0.29955511
## AAACCCAAGTCATCCA-1 Microglia -4.1035505 -0.76945994   -0.36419246
## AAACCCAAGTCGGCCT-1 Microglia 10.8851187  2.14596437    0.82758593
## AAACCCACAAGAGGTC-1 Microglia  2.9593327  0.05956462    0.29261936
##                    DAM.sig_top100 DAM.sig_top250 DAM.sig_top500    m1.blue
## AAACCCAAGAAGTCCG-1    -0.05322557   0.0006620787      0.0917550 0.10662676
## AAACCCAAGAGAAGGT-1     0.01828038   0.1322962948      0.2170636 0.11369759
## AAACCCAAGCACTCAT-1     0.28877250   0.2463909542      0.3972866 0.08340758
## AAACCCAAGTCATCCA-1    -0.24372908  -0.0888669115      0.1221884 0.15953562
## AAACCCAAGTCGGCCT-1     0.58889142   0.5116411395      0.5554119 0.06910936
## AAACCCACAAGAGGTC-1     0.22824103   0.2793449722      0.3333028 0.04203820
##                     m2.purple m3.darkorange   m4.yellow m5.yellowgreen
## AAACCCAAGAAGTCCG-1 0.19167658   0.134384240 -0.06013164    -0.01712984
## AAACCCAAGAGAAGGT-1 0.17418624   0.113741642  0.10167319     0.02485532
## AAACCCAAGCACTCAT-1 0.11077992   0.197485327  0.21626831    -0.05288063
## AAACCCAAGTCATCCA-1 0.17844982   0.003841682 -0.03585387    -0.08453258
## AAACCCAAGTCGGCCT-1 0.07126132   0.240925300  0.08158639     0.02887132
## AAACCCACAAGAGGTC-1 0.11175430   0.242583580  0.07825132     0.18873265
##                     m6.green m7.darkgreen m8.turquoise
## AAACCCAAGAAGTCCG-1 0.3442975 -0.038220365  -0.03207072
## AAACCCAAGAGAAGGT-1 0.5728791 -0.031307089  -0.06703069
## AAACCCAAGCACTCAT-1 0.8033512 -0.016257571   0.07148222
## AAACCCAAGTCATCCA-1 0.4835077 -0.054728716  -0.01170676
## AAACCCAAGTCGGCCT-1 0.9344805  0.089015706   0.08991283
## AAACCCACAAGAGGTC-1 0.8320493  0.005748232  -0.03052330
#saveRDS(GEX.seur,"I:/Shared_win/projects/202310_Spp1DAM/forGEO/seur_obj/Spp1GFP_Modelx3.final.seur_obj.rds")